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SECTION  I 


INTRODUCTION 

In  this  report,  we  shall  study  the  response  of  a jet-engine  fan  blade  due 
to  impacts  by  foreign  objects,  such  as  birds,  ice  and  runway  stones.  The  fan 
blade  will  be  treated  as  a laminated  composite  plate.  Our  goal  is  to  calculate 
or  estimate  the  stress  level  and  damage  threshold  in  the  blade  and  present  the 
resulting  data  in  a manner  that  can  be  of  use  to  the  designer. 

The  report  is  divided  into  two  parts:  the  first  involves  structural  re- 
sponse, and  presents  a procedure  for  constructing  design  curves.  The  second 
part  concerns  local  response  and  the  development  of  an  anisotropic  finite-differ- 
ence code.  In  the  structural-response  part  of  the  report,  we  shall  first  review 
a recently  developed  approach  of  treating  the  parameters  in  an  impact  problem 
in  a particular  non-dimensional  form.  A single  design  curve  can  then  be  applied 
to  the  same  type  of  structure  under  various  impact  conditions.  The  previously 
published  paper  covered  only  simply  supported  beams,  but  in  this  report  the 
approach  will  be  extended  to  simply  supported  plates,  cantilever  plates  and  the 
effect  of  an  impact  by  a soft  object,  such  as  a bird  or  water.  A single  design 
curve  is  presented  for  each  structure  and  the  procedure  for  generating  additional 
design  curves  by  either  a numerical  or  experimental  method  is  also  discussed. 

For  the  local-response  calculation,  either  finite-element  or  finite-differ- 
ence methods  can  be  used.  The  finite-element  method  is  more  convenient  for 
treating  geometrically  irregular  structures;  on  the  other  hand,  the  finite-dif- 
ference methods  have  a long  history  of  solving  wave-propagation  problems  and 
there  are  many  codes  available  for  adaption  to  our  particular  problem.  Our  com- 
posite fan  blade  can  be  smeared  into  an  equivalent  homogeneous  anisotropic 
material.  The  geometry  of  the  structure  is  then  greatly  simplified  and  there- 
fore the  problem  can  be  handled  conveniently  by  the  finite-difference  method. 


There  are  many  existing  finite-difference  computer  codes,  for  example, 

CEL  [1],  TENSOR  [2],  HELP  [3]  and  HEMP  [4].  All  of  these  codes  are  restricted 
to  isotropic  materials.  We  have  modified  the  HEMP  code  to  include  anisotropic 
elastic  constitutive  relations  and  resulted  with  the  ANEL  code  (Anisotropic 
Elastic) . A few  sample  impact  problems  were  solved  by  ANEL  and  the  results 
are  compared  with  experimental  data  and  other  existing  solutions. 

In  Section  II  the  Design  Curve  approach  is  presented;  in  Section  III  the 
ANEL  code  is  described.  These  two  sections  are  presented  independently. 


1.  Non-dimensional  Parameters  for  Structural  Response 
In  Ref.  [5],  the  method  for  constructing  a design  curve  for  predicting 
the  response  of  a given  structure  to  impact  is  presented.  This  curve,  which 
is  in  terms  of  a generalized  strain  e,  and  the  mass  ratio  M (impactor  mass/ 
structure  mass)  can  be  applied  to  all  structures  of  the  same  type  but  of 
different  dimensions  and  various  materials.  It  also  covers  impactors  of  differ- 
ent mass  and  velocity.  Simply  supported  beams  were  used  as  an  example  in  [5], 
although  the  same  approach  can  be  applied  to  other  types  of  structures. 

In  demonstrating  the  method.  Ref.  [5]  investigated  six  analytical  models. 


Three  of  these  were  lumped-parameter  models  consisting  of  discrete  masses  and 
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springs;  the  other  three  were  distributed-mass  models,  which  took  the  beam 
vibration  modes  into  consideration.  Some  of  these  models  neglected  the  contact 
force  and  treated  the  beam  with  the  attached  impactor  mass  as  a free  transient- 
vibration  problem.  Others  included  the  contact-force  term,  traced  the  impactor 
motion  distinctly,  and  treated  the  beam  as  a forced-vibration  problem.  It  was 
shown  that  for  the  models  neglecting  the  contact  force,  the  generalized  strain 
e is  a function  of  mass  ratio  M only  and  is  independent  of  other  parameters. 

For  the  two  models  which  included  the  contact  force,  e depends  on  M and  other 
parameters,  but  the  dependence  on  these  other  parameters  is  weak,  and  a single 


e vs.  M curve  can  still  give  an  approximate  representation  of  all  impact 


situations . 

In  this  report,  we  shall  adopt  two  of  these  six  models;  the  Timoshenko 
solution  (or  model)  and  a single-degree-of-f reedom  model.  The  Timoshenko  model 
applied  to  simply  supported  plates  will  be  reviewed  in  detail  in  a later  section. 
Here,  we  shall  outline  the  single-degree-of-freedom  model. 
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In  this  model,  the  impactor  and  the  structure  are  considered  attached 
together  immediately  on  contact  as  a lumped  mass  m,  and  the  deflection  of 
the  structure  is  governed  by  an  equivalent  spring  of  spring  constant  K^. 

The  initial  velocity  of  the  combined  mass,  v^,  may  be  determined  by  two 
methods:  one  is  based  on  the  conservation  of  momentum,  the  other  on  the  con- 
servation of  energy.  In  the  first  case 

m = + em^  (1) 

where  em^  is  an  "equivalent"  mass  of  the  structure.  The  equivalent  mass  may 
be  obtained  by  matching  the  kinetic  energy  of  the  equivalent  mass  traveling 
at  the  velocity  of  the  impact  point  in  the  structure  and  the  total  kinetic 
energy  of  the  structure,  assuming  that  the  deflection  mode  shape  is  the  static 
deflection  curve.  For  simply  supported  beams,  e = 17/35.  By  equating  the 
momentum  of  the  impactor  before  impact,  m2  v,  with  the  momentum  of  the  combined 
mass  m,  we  have  (M  = 

Vq  = v/  (1  + eM)  . (2) 

Note  that  according  to  Eq.  (2)  the  energy  before  and  after  impact  is  not  con- 

2 2 
served,  or  mv^  < n^v  . 

In  the  energy-conserved  case,  we  assume  e = 0,  and 


vQ  = v 

As  a result,  the  kinetic  energy  in  the  impactor  is  conserved,  and  will  be 


entirely  converted  to  strain  energy  in  the  structure.  For  small  values  of  M, 
the  difference  between  these  two  approaches  is  small. 

The  spring  constant  is  the  force  per  unit  deflection,  with  respect 
to  a force  acting  at  the  impact  point  in  the  direction  of  the  impactor 


velocity,  or 


P " K1  W1 
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After  acquiring  the  initial  velocity  Vq,  the  mass-spring  system  is 
assumed  to  perform  a free  vibration.  The  maximum  deflection  is  then 


W1  = vo  ‘ 

max  * 1 


In  replacing  the  structure  by  a mass-spring  system,  it  is  implied  that  only 
the  first  mode  of  the  structure  is  retained  and  that  its  mode  shape  is  the 
same  as  the  static  deflection  distribution  under  a concentrated  force. 

Next,  we  shall  consider  the  most  critical  strain  in  the  structure 
which  occurs  at  a known  point  2.  Assume  that  a relation  between  z^  and  w^ 


can  be  found, 

» 


e,  = — wi 

2 d12  1 


where  d ^ can  be  determined  from  the  static  deflection  distribution,  w^(x), 
or  from  static  measurement.  Combining  eq.(5)  and  eq.  (6),  we  obtain  the  maxi- 


mum strain  as 


'2  max  d 


or,  combining  Eq.  (1),  (2)  and  (7), 


'o  /jrT 

.2  ih 


2 max  d 2 |f  (l+eM)K1 

The  next  step  of  defining  a generalized  strain  e is  best  illustrated  by 
considering  a specific  type  of  structure.  For  a simply  supported  beam  impacted 


at  middle  span,  we  have  [5] 


d 2 = L /6h 


K = 48EI/L  . 
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Therefore, 

e 

max 

where  Cq  = /E/p  , the  velocity  of  longitudinal  waves  in  a bar.  Defining  the 
generalized  strain  as 


v h FT" 
= co  k V4 


M ( 1+eM) 


(ID 


we  obtain  finally, 


e 


e 


max 


fo  k 

v h 


(12) 


e 


1 

M(l+eM) 


(13) 


which  is  the  equation  of  the  previously  mentioned  design  curve. 

The  corresponding  expressions  for  d^»  and  e for  simply  supported 
orthotropic  plates,  cantilever  beam,  and  cantilever  plates  will  be  given  in 
later  sections. 

2.  Numerical  Solution  of  Simply  Supported  Plate 

The  engine  fan  blade  will  be  simulated  by  a cantilever,  plate  later  in 
this  report.  Since  it  is  not  convenient  to  express  the  deflection  of  a 
cantilever  plate  in  explicit  form,  we  shall  depend  on  the  one-degree-of- 
freedom  model. 

In  order  to  study  the  approximation  involved  in  using  the  one-degree-of- 
freedom  model,  we  shall  first  analyze  a simply  supported  plate,  which  possesses 
a more  exact  solution,  the  Timoshenko  solution.  The  simply  supported  plates 
will  also  be  used  to  study  the  approximation  involved  in  simulating  a high- 
speed bird  impact  by  a low-speed  hard-object  impact.  Both  the  simple  one- 
degree-of-f reedom  model  and  the  Timoshenko  model  will  be  used, 
a.  One-degree-of-f reedom  Model 

A simply  supported  rectangular  orthotropic  plate  under  a concentrated 
load  at  the  center  is  considered.  The  maximum  deflection  occurs  at  the  center 
of  the  plate  and  the  maximum  strain  on  the  surface  opposite  from  the  applied 
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load.  It  can  be  shown  (see  Appendix  A)  thac  under  a mild  assumption,  the  spring 


constant  of  the  plate  is 


K,  = — 
1 w 


where 


A 

it  Dx 

4 ab(§)2  fjW 


co  ou 

Mn)  - l l 

m=l,3,5  n=l,3,5  mn 


= fzVJ^L 


4,22^42  n7s 

C = m + 2m  n n + n n v1'-* 

mn 

The  strain-deflection  constant  dl2  relating  the  strain  in  the  x-direction  with 
w^  is 

e . 7T2h  f.(n) 

_*  = _i_  = l (18) 

W1  d12  2 ab(J)  fx(n) 

where 


Mn)  l l r~ 

m=l,3,5  n=l,3,5  mn 

Combining  Eq.  (14),  (18)  with  Eq.  (8),  and  assuming  e = 0,  we  have 


e.  = f0(h)  h v 
"max  2 


’ P h T 

. Dx  M f^(n)  J 


If  we  define  the  generalized  strain  as 


E = £ 

x x max 


where  c^  is  the  speed  of  flexural  waves  in  the  x-direction -^D^TT  , then 


'x  r u £ t t i 2 
[M  f1(n) J 


Similarly, 


f2(n_1) 


-U.l/2 


y [m  fx(n  )] 
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The  vs.  M curve  based  on  eq-  (22)  for  a particular  plate  will  be  given 
later. 

b.  Timoshenko  Solution  of  Transverse  Plate  Impact 
Timoshenko's  approach  for  solving  transverse  impact  problems  on 
beams  [6]  by  coupling  Hertz's  law  of  contact  with  the  Euler  beam  equation  has 
been  extended  to  the  case  of  a simply  supported  isotropic  rectangular  plate 
by  Karas  [7]  and  more  recently  to  the  case  of  a simply  supported  anisotropic 
laminated  plate  by  Sun  and  Chattopadhyay  [8].  In  this  section  we  review  the 
salient  points  of  the  latter  solution  and  show  how  it  may  be  used  to  develop 
an  impact  design  curve  for  anisotropic  plates. 

Solution  Method  * 

The  anisotropic  plate  equations  of  Whitney  and  Pagano  [9]  are  used  to 
predict  the  motion  of  the  plate.  Neglecting  the  effects  of  rotatory  inertia, 
the  deflection  of  a symmetric  cross-ply  laminated  plate  due  to  a centrally 


applied  force  F(t)  may  be  expressed  as 


w(x,y , t)  = 


l l ^ 


nrhi-2 
2 rrt 


m,  “ “ w 

1 m n mn 


. irnrx  . niry 

•sin  sm 

a b 


F(x)sin  ^(t-xjdxj 


m.n  = 1,3,5  . . °° 


where  m, , a,  and  b are  the  mass  and  planar  dimensions  of  the  plate,  and  w 

1 mn 

are  natural  frequencies  dependent  on  m and  n and  the  properties  of  the  plate. 
Hertz's  law  of  contact  is  assumed  to  hold 

F = k2  a3/2  (25) 

where  k2  is  a constant  and  a is  the  indentation  of  the  impactor  relative  to 


the  plate  surface,  or 


a = w2  - w 


where  w2  is  the  impactor  displacement. 


Newton's  law  applied  to  the  impactor  is 


w2  = vt 


- — r r 

m2  J 0 0 


F dt  dt 


These  four  equations  may  be  combined  into  a single  nonlinear  integral 
equation  in  terms  of  the  contact  force  F between  the  plate  and  the  impactor, 


til 


= vt  - 


2 •>  0 •'0 


F dt  dt 


00  OU  f £ 

— y y — — F(x)sin  to  (t-x)dx 

m.  .L _ c , (o  Ja  11111 

T.  m=l,3,5  n=l,3,5  mn  J0 

which  is  solved  numerically  by  applying  the  small- increment  method  suggested 
by  Timoshenko,  in  which  the  contact  force  F is  assumed  to  be  constant  during 
any  time  increment  Ax.  Expanding  the  above  integrals  to  calculate  the  force 
F^  during  the  itk  time  interval,  we  obtain 


F.  -I  2/3 


ry 

. k? . 


= vnAx  - 


where 


^ fj  - 1 k Fl  '*** 

s =y  y — 

i-i+1  L L 2 


-j+l  L L 2 

J m n a) 

mn 


•{cos[to  (i-j)Ax]  - cos[oj  (i-j+l)Ax]},  m,n  = 1,3, 5,..., 00  (30) 

mn  mn 

If  the  contact  force  is  approximated  as  a linear  continuous  function  of 

time,  with  an  average  value  of  F^  during  the  ith  time  step,  then 
i 

l Fj  = 2[(n-l)F7  + (n-2)(F2  - F^ 

+ (n-3)(F3  - F2  + Fj)  +•••+  (n-jHFj  - F x + Fj-2 


•+  F. )+•••+  (F  . - F + F 
— 1 n-1  n-2  n-3 


•±  VI  + 3(Fn-  Fn-l  + Fn-2 


•+  Fx) 
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For  computing  the  solution  of  Eq.(29),  Sun  and  Chattopadhyay  suggested 
a recursion  method,  but  we  have  found  that  such  a time-saving  approximation  is 
not  necessary.  A listing  of  our  FORTRAN  computer  program  is  presented  in 
Appendix  B. 


Construction  of  Design  Curve 

As  can  be  seen  from  the  nonlinear  nature  of  Eq. (28),  the  structural 
response  of  the  plate  to  impact  is  dependent  on  more  parameters,  not  only  the 
mass  ratio  M = m^/m^.  However,  we  will  demonstrate  that,  for  plates  of  given 
aspect  ratio  (a/b)  and  given  anisotropy  ratio  (°22^11^ ’ dependence  °f  the 


plate  response  on  these  other  parameters  is  weak,  and  that  one  family  of  e vs. 

M curves  gives  a good  approximation  of  all  impact  cases. 

To  determine  the  effect  of  these  quantities  on  impact  response,  a para- 
metric calculation  was  carried  out,  with  different  values  of  the  following 
parameters:  impact  velocity  v,  contact  stiffness  k2»  plate  bending  stiffness 
matrix  [D^],  and  mass  ratio  M.  The  actual  magnitudes  of  these  quantities  used 
and  the  calculated  values  of  the  generalized  strains  ex  and  are  all  summar- 


ized in  Table  I. 


An  inspection  of  this  table  indicates  that  for  any  given 


value  of  M,  the  maximum  difference  between  any  two  values  of  e , or  of  e , is 
only  12%. 

This  solution  may  be  used  to  construct  an  impact  design  curve  by  plotting 
points  on  an  e vs.  M graph  (Figure  1).  These  points  alone  can  be  connected 
by  a single  curve  which  may  be  used  as  a design  curve.  Note  that  such  a curve 
would  not  be  very  different  from  the  curve  corresponding  to  the  one-degree-of- 
freedom  model  discussed  in  an  earlier  section. 
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Table  1 PARAMETRIC  STUDY  OF  SIMPLY- SUPPORTED  PLATES  BY 
TIMOSHENKO  SOLUTION,  (v  = 2.45  m/sec,  k2  = 1.744xl019  N/m3/2  ) 


\ 


i=S5ss=::. 


iipiHfilllill! 


:«:»:! 


Figure  1.  Central  Impact  of  Simply  Supported  Rectangular  Orthotropic  Plate 


c.  Bird  Impact 

Under  high-speed  impact,  the  body  of  a bird  responds  as  a fluid  with 
negligible  strength.  In  this  section,  we  shall  treat  the  bird  as  a fluid 
and  compute  the  stress  induced  in  the  blade  due  to  fluid  impact.  The 
Timoshenko  solution  for  the  response  of  a simply  supported  plate  to  a solid 
impactor  was  modified  to  treat  the  loading  of  a fluid  impact.  Instead  of  a 
Hertzian  contact  force,  the  force  F exerted  on  the  plate  is  given  by 

F = (32) 

where  $ is  the  density  of  the  impacting  fluid,  A the  cross-sectional  area 

of  the  fluid  column,  and  v is  the  relative  velocity  between  the  fluid  column 

R 

and  the  plate  impact  point.  It  is  assumed  that  the  fluid  column  cannot  sus- 
tain axial  compression,  and  the  fluid  velocity  v remains  constant  until  the 
complete  fluid  column  is  consumed.  The  relative  velocity  vR  is  then 

vR  = v - w (j,  y)  (33) 

Substituting  Eqs.  (32)  and  (33)  into  (24)  we  obtain  the  governing 
equations  for  w;  this  equation  is  then  solved  numerically. 

The  contact  duration  T,  which  is  the  time  required  to  completely  consume 
the  column,  can  be  approximated  as 

T = y (34) 

where  L is  the  length  of  the  fluid  column.  For  t > T,  the  force  F is  zero 
and  the  plate  will  continue  free  vibration. 

Calculations  were  made  for  two  velocity  regimes;  the  low-speed  regime 
was  around  6 m/sec,  and  the  high-speed  regime  larger  than  60  m/sec.  The  low- 


speed  impact  cases  were  studied  in  order  to  compare  with  the  results  of  the 
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gelatin  impact  experiments.  The  contact  time  is  long  for  low-speed  impact, 

and  the  fluid  loading  force  fluctuates  a very  small  amount  from  that  given 
2 2 

by  i)>Av  , instead  of  4>Av  . In  other  words,  the  motion  of  the  plate  does  not 

K 

affect  the  contact  force  too  much,  as  shown  in  Figure  2.  The  strain  is 

approximately  the  same  as  that  produced  by  a static  force  of  the  magnitude 
2 

<t>Av  , Figure  3.  The  results  will  be  discussed  further  in  the  section  on 
experimental  work. 

The  high  speed  impact  cases  were  studied  in  order  to  simulate  the  true 
high-speed  bird  impact.  The  contact  time  is  short  and  the  plate  velocity,  w 
becomes  large  so  that  its  motion  influences  the  contact  force.  A parametric 
study  was  undertaken  in  which  M,  the  ratio  of  the  plate  mass  to  impactor  mass, 
A,  the  cross-sectional  area  of  the  fluid  column,  and  v,  the  impactor  velocity, 
were  varied.  Results  of  this  study  are  presented  in  Table  2 and  Figure  4. 

As  can  be  seen,  points  plotted  in  coordinates  of  generalized  strain  e vs.  mass 
ratio  M are  grouped  in  a single  band.  For  any  value  of  M,  the  largest  spread 
of  values  of  e is  only  34%,  which  suggests  that  a curve  through  these  points 
may  be  useful  as  a design  tool  for  predicting  the  response  of  other  simply 
supported  plates  to  fluid  impact.  Further,  by  conducting  a few  impact  tests 
on  actual  fan-blade  structures,  a similar  curve  can  be  constructed  to  serve 
as  a design  guide  for  a wide  range  of  impact  cases. 


train  vs.  Time  Low  Spee< 


Table  2.  PARAMETRIC  STUDY  OF  HIGH-SPEED  FLUID  IMPACT  MODEL 


M 

v,  m/sec 

A 2 
A,  mm 

e mm/ mm 

e 

1.92 

60.96 

5806. 

0.031 

0.56 

1.49 

60.96 

5806. 

0.032 

0.58 

91.44 

5806. 

0.058 

0.71 

121.92 

5806. 

0.054 

0.49 

152.40 

5806. 

0.074 

0.54 

274.32 

5806. 

0.140 

0.56 

0.96 

137.16 

5806. 

0.121 

0.97 

0.75 

91.44 

5806. 

0.065 

0.79 

121.92 

5806. 

0.105 

0.95 

152.40 

5806. 

0.137 

0.99 

0.64 

137.16 

7903 

0.133 

1.07 

0.49 

91.44 

7903 

0.083 

1.01 

121.92 

7903 

0.125 

1.14 

152.40 

7903 

0.171 

1.24 

60.96 

17419 

0.066 

1.19 

152.40 


17419 


0.131 


0.95 


; 

:::an 


3.  Experimental  Approach 

Several  series  of  impact  experiments  were  conducted.  All  the  experi- 
ments used  composite  plates  of  the  same  material  and  layup,  but  with  different 
supports.  One  series  involves  a simply  supported  plate,  the  other  cantilever 
plates.  The  simply  supported  plate  was  impacted  transversely  both  by  steel  im- 
pactors  and  by  gelatin.  The  cantilever  plates  were  impacted  by  steel  and 
aluminum  impactors  in  the  transverse  direction  and  on  the  edge  at  various  angles 
of  incidence. 

a.  Specimens 

All  plate  specimens  were  fabricated  from  Hercules,  Inc.,  AS3501 
graphite-epoxy  prepreg  tape,  at  a layup  of  [ (0/45/0/-45)g/0]s.  The  dimensions 
and  mass  of  the  simply  supported  plate  were  349x171x7.1  mm  (a  x b x h)  and 
0.660  kg.  Two  cantilever  plates  were  tested:  (1)  330x178x7.1  mm,  0.647  kg., 
and  (2)  330x89x7.1  mm,  0.324  kg.  Plate  (2)  was  cut  from  plate  (1).  These 
dimensions  and  masses  include  only  the  supported  areas  of  the  plates;  an  additional 
4-mm  margin  protruded  beyond  the  supports  around  the  simply  supported  plate,  and 
24  mm  of  the  cantilever  plates  extended  into  the  clamping  device. 

The  simply  supported  plate  and  the  cantilever  plate  were  transversely 
impacted  by  blunt  (25. 4-mm  contact  radius)  steel  cylindrical  impactors  using 
a drop-test  machine.  The  cantilever  plate  was  also  tested  with  edge  impacts. 

b . Simply  Supported  Plates 

The  simply  supported  plate  was  rested  upon  a rectangular  frame  of  round- 
edged  blades  and  impacted  at  the  center  by  various  impact  masses  at  several 
velocities.  The  strain  in  the  plate  directly  opposite  the  impact  point  was 
measured  using  a Micro-Measurements,  Inc.,  type  EA-41-125AD-120  strain  gage. 


and  recorded  using  an  Ellis  BAM-1  bridge  amplifier  and  Tektronix  Type  565 


oscilloscope.  Impact  velocity  was  determined  from  the  drop  height.  Data 
from  these  experiments  by  steel  impactor  are  summarized  in  Table  3.  This 
result  is  also  included  in  Figure  1. 

c.  Simulation  of  Bird  Impact  by  Gelatin 

An  attempt  was  made  to  simulate  high-speed  bird  impacts  by  low-speed 
impacts  produced  by  dropping  a gelatin  mixture  onto  the  simply  supported  plate. 
The  strain  in  the  plate  was  measured  in  the  same  manner  as  described  above  for 
the  hard-object  impact  experiments.  The  gelatin  mixture  was  prepared  from 
food-type  gelatin  (e.g.,  Jello)  according  to  package  directions  but  with  75% 
extra  water  added. 

After  much  experimentation  with  various  devices  for  dropping  the  gelatin 
onto  the  plate,  fairly  repeatable  results  were  finally  obtained  by  the  following 
procedure.  The  gelatin  mixture  was  prepared  (set)  inside  lengths  of  2.82-inch 
(I.D.)  plastic  pipe.  A smooth  flat  plate  was  firmly  pressed  against  the  bottom 
end  of  the  pipe  to  retain  the  gelatin.  The  pipe  was  held  above  the  plate 
specimen  and  the  retaining  plate  was  suddenly  pulled  aside,  allowing  the  gelatin 
to  fall.  A sheet  of  aluminum  foil,  situated  a small  distance  above  the  plate 
specimen,  made  contact  with  the  plate  on  impact  of  the  gelatin,  completing  an 
electrical  circuit  and  triggering  the  oscilloscope  for  strain  measurement. 

In  the  course  of  these  experiments,  it  was  observed  that  the  impact  dura- 
tion was  longer  than  the  calculated  period  of  the  lowest  mode  of  plate  vibration. 
However,  in  an  actual  high-speed  bird  impact,  the  impact  duration  is  character- 
istically much  shorter  than  the  fundamental  period.  This  large  discrepancy  in 
the  kinematics  of  the  two  cases  indicates  that  low-speed  fluid  impact  experi- 
ments cannot  be  used  to  even  approximately  simulate  high-speed  bird  impact. 
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d.  Cantilever  Plates 

For  all  impacts  on  the  cantilever  plate,  the  specimen  plates  were 
held  in  a jig  that  provided  25  mm  of  clamping  depth.  For  the  edge  impact 
experiments,  the  jig  could  be  rotated  to  various  angles  of  incidence. 

Two  series  of  experiments  were  conducted  for  the  cantilever  plates, 
transverse  impact  and  edge  impact.  These  will  be  described  separately  below. 

(1)  Transverse  Impact 

Cantilever  plate  (1)  was  impacted  transversely  in  a manner  similar  to  the 
simply  supported  plate,  at  a point  76  mm  from  the  free  end.  Bending  strain  in 
the  plate  was  measured  at  a point  25.4  mm  from  the  clamped  support  using  a 
strain  gage  (Figure  5).  The  results  of  these  tests  are  presented  in  Table 
4 . 

Cantilever  plate  (2)  was  also  failure-tested  by  gradually  increasing  both 
the  impact  mass  and  drop  height  until  failure  was  visually  detected.  During 
each  test,  the  strain  in  the  plate  was  measured  by  a strain  gage  mounted  25  mm 
from  the  clamped  end.  The  plate  finally  failed  when  impacted  by  a 6.08-kg 
lead  block  at  5.99  m/sec.  Failure  occurred  in  a zone  about  3 to  5 cm  away 
from  the  clamped  end,  and  apparently  resulted  from  the  high  bending  stresses 
produced  in  this  area.  Strain-measurement  results  are  also  presented  in 
Table  4. 

(2)  Edge  Impact 

This  experimental  study  was  the  first  phase  of  an  experimental  analytical 
study  which  had  as  its  goal  the  development  of  an  analytical  tool  to  predict 
the  transient  strain  response  of  a turbine-blade-like  structure  subjected  to 
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Specimen  failed  during  this  test 


impact  on  its  edge  at  various  angles  of  incidence.  However,  due  to  a change 


in  scope  of  this  grant  the  analytical  phase  of  this  study  was  never  initiated. 

The  experimental  results  are  presented  here  for  reference  purposes  in  case 
this  proposed  study  is  reactivated  in  the  future. 

This  series  of  impact  tests  involved  impacting  the  cantilever  specimen 
plate  (previously  described)  on  its  edge  by  a 0.85"  (21.6mm)  diameter  by 
1.0"  (25.4mm)  long  aluminum  cylindrical  impactor.  The  impactor  was  accelerated 
by  a low-pressure  air  gun  to  a velocity  of  40  fps  (12.2  m/sec);  the  point  of 
impact  on  the  specimen  plate  was  at  its  80%  length.  During  the  tests  the 
cantilever  plate  was  rotated  such  that  the  incident  angles  of  impact  included 
in  this  study  were  0°,  22.5°,  45°,  67.5°  and  90°.  Figure  6 is  a photo- 

graph of  the  experimental  set-up  used  in  these  tests.  Figure  7 identifies 
the  locations  of  the  strain  gages  and  point  of  impact  of  the  cylindrical  impactor. 

At  least  three  sets  of  reproducible  strain  traces  were  recorded  for  each 

*• 

angle  of  incidence  impact.  Figures  8,  9,  and  10  are  typical  strain 
histories  of  gages  4 and  9;  for  the  0°  impact  these  traces  are  identical 
and  for  the  90°  impact  they  are  identically  out  of  phase.  Table  5 in- 
cludes the  values  of  the  average  of  the  peak  strains  at  each  gage  for  each 
angle  of  incidence.  Also  included  in  this  Table  is  the  experimental  value  of 
the  wave  velocity  based  on  strain  initiation  at  gages  1 and  4 as  the  wave 
propagates  into  the  plate.  As  anticipated,  the  value  of  the  wave  velocity 
decreases  as  the  impact  changes  from  an  inplane  mode  (0°)  to  a transverse  mode 
(90°).  The  theoretical  values  of  the  inplane-  and  shear-wave  velocities  for 
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Note: 


Gage  9 is  on  top  of  plate;  all  others  on  bottom. 
All  dimensions  in  millimeters. 


Figure  7 . 


Specimen  Geometry  and  Strain  Gage  Locations 


gage  1 
gage  2 

gage  3 
gage  4 
gage  9 


gage  1 

gage  5 
gage  6 
gage  7 


Horizontal  Scales:  50  usee  per  major  division 
Vertical  Scales:  0.1%  strain  per  major  division 

Figure  8 . Typical  Oscilloscope  Traces  for  Edge- 

Impact  Experiments  at  0°  Angle  of  Incidence. 


\ 


28 


Horizontal  Scales:  50  nsec  per  major  division 
Vertical  Scales:  0.1  % strain  per  major  division 


Typical  Oscilloscope  Traces  for  Edge- 
Experiments  at  45°  Angle  of  Incidence 


j 


gage  1 
gage  2 

gage  3 
gage  4 

gage  9 


gage  1 

gage  5 
gage  6 
gage  7 


Horizontal  Scales:  50  ysec  per  major  division 
Vertical  Scales:  0.1%  strain  per  major  division 

Figure  10  . Typical  Oscilloscope  Traces  for  Edge-Impact 

Experiments  at  90°  Angle  of  Incidence. 
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4.  Design  Curves  for  Cantilever  Plate 

Since  the  exact  solution  of  the  static  deflection  of  orthotropic  rectan- 
gular cantilever  plates  is  not  known,  we  cannot  obtain  explicit  analytical 
expressions  for  and  d^*  Instead,  we  shall  handle  this  in  two  approximate 
ways.  The  first  is  to  treat  the  plate  as  a cantilever  beam,  and  use  the  de- 
flection function  of  the  beam  to  obtain  and  d^*  The  other  method  is  to 
measure  K.^  and  d^  directly  from  static  tests. 

We  shall  consider  the  plate  as  a cantilever  beam  of  length  L under  a 
concentrated  load  P applied  at  a point  x = d as  shown  in  Figure  11.  The 
expression  for  the  strain  at  a point  x = c will  be  derived.  The  deflection 
distribution  of  a cantilever  beam  subjected  to  such  a load  is 


wx(x)  = 


(3d-x) 


( 3x-d ) 


0 < x < d 


d < x < L 


At  x = d. 


so  that 


,,,  Pd 
Wl(d)  = 2EI 


The  bending  strain  in  the  beam  at  x = c may  also  be  determined  from  the 
static  deflection  curve. 

2 

ll  3 w 

E2  " 2 , 2 
3x 


(c  < d) 
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Therefore 


12  " h(d-c) 

Following  the  definition  of  e from  Eq.(12),  we  have 

1/2 


(39) 


f1  " d)  ^4M(d/L)_ 


(40) 


In  this  case  e is  a function  not  only  of  M but  also  of  the  parameters 
(d/L)  and  (c/d).  Note  that  e is  a maximum  at  c = 0,  i.e.,  at  the  supported 
end . 


Another  way  of  determining  and  d^  for  the  cantilever  plate  is 

direct  static  measurement  of  P,  w^ , and  £2*  Referring  to  Figure  11,  if 
a load  P is  applied  at  x = d,  and  the  deflection  w^  at  x = d and  the  strain 
£2  at  x = c are  measured,  then  the  values  of  and  d^  can  be  calculated 
from  Eqs.  (4)  and  (6). 

We  have  made  static  measurements  on  the  specimen  shown  in  Figure  5. 
The  results  are  presented  in  Table  6 . The  experimental  values  of 

and  d^2  ate 

K = 49700  N/m 

d = 12.6  m 


Applicat/on  of  Design  Curves  to  Turbine  Blades 

The  design  curve  for  cantilever  plates  and  beams  can  be  used  to  deter- 
mine the  resistance  of  a turbine-engine  fan  blade  to  structural  failure 
caused  by  a particular  impact.  This  will  be  demoostrrated  by  means  of  an 
example. 

For  a blade  similar  to  the ' composite  cantilever  plate  specimen  de- 
scribed in  section  11-3'  (330  x 178  x 7.1  mm,  0.647  kg),  it  is  desired  to  find 


f 


Table  6.  STATIC  TESTS  ON  CANTILEVER  COMPOSITE  PLATE  SPECIMEN 


Force 


Deflection 


Strain 


lb 

N 

.001  in 

ram 

ys 

0.0 

0.0 

0.0 

0.000 

0 

5.0 

22.2 

15.6 

0.396 

36 

10.0 

44.5 

33.0 

0.838 

72 

15.0 

66.7 

48.7 

1.237 

106 

20.0 

89.0 

65.0 

1.651 

140 

23.1 

102.6 

78.0 

1.981 

166 

25.1 

111.5 

84.5 

2.146 

178 

30.1 

133.7 

101.4 

2.576 

216 

35.1 

156.0 

119.2 

3.028 

252 

40.1 

178.2 

137.0 

3.480 

288 

44.9 

199.6 

153.2 

3.891 

316 

45.1 

200.4 

154.3 

3.919 

322 

49.9 

221.9 

170.8 

4.338 

352 

54.9 

244.1 

189.3 

4.808 

388 

59.9 

266.3 

208.6 

5.298 

428 

61.9 

275.2 

216.1 

5.489 

442 
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the  structural-failure  impact  velocity  for  a 0.54-kg  bird  striking  the  blade 
76  mm  from  the  tip.  From  static  flexure  tests,  it  is  determined  that  the 
flexural  rigidity  El  is  408  N-m  (therefore  a = /EI/pA  = 14.4  m /sec),  and 
that  the  blade  fails  in  bending  near  the  supported  end  when  the  outer-fiber 
tensile  strain  reaches  the  value  of  e = 0.008.  With  these  data  the  design 
curve  for  the  blade  may  be  constructed.  Treating  the  blade  approximately  as 
a cantilever  beam,  the  equation  for  the  curve,  Eq.(40),  may  be  written  as 


where 


=i-  M)VS 

■ G ■ 254)  ■/ 


3(330) 

4(254) 


= 0.866 

For  a mass  ratio  M = 1.2,  the  generalized  strain  is  e = 0.79.  Therefore,  the 
normal  impact  velocity  which  will  just  initiate  failure  is 


(0.008) (14.4) 
(0.0071) (0. 79) 

= 20.5  m/sec 


A design  curve  for  the  blade  may  also  be  constructed  by  means  of  static 
strain  and  deflection  measurements.  In  this  case,  the  constants  and  d^ 
are  determined  experimentally,  as  discussed  in  the  previous  section.  The 


values  are  found  to  be 


K.  = 49700  N/m 

d^  = 12.6  m 

for  a static  force  applied  at  the  impact  point  and  for  strain  measured  near 


I 

I 

* 


the  supported  end.  Then  by  Eq.(8),  using  e = 0 (conservation-of-energy 

model)  we  have 


= (12.6)  (0.008) 
= 30.6  m/sec 


The  equation  for  the  corresponding  design  curve  is 

- 

£ vfi 


where 


=1 IV 


d12h 


0.582 


This  equation  may  also  be  used  to  determine  the  critical  impact  velocity  as 
described  above.  These  two  equations  for  the  design  curve  are  plotted  in 
Figure  12. 

Design  Curve  in  Terms  of  Dimensionless  Velocity 

The  impact  design  curve  may  be  presented  in  another  useful  form.  A 
dimensionless  velocity  is  defined  as 


v = 


v 


h 


2 


where  is  the  value  of  strain  for  which  the  structure  fails.  Note  that  when 
v has  the  value  of  the  critical  or  failure  velocity  vc>  the  dimensionless 
velocity  is  the  reciprocal  of  the  generalized  strain  e. 


v 

c 


t 
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Figure  12.  Impact  Design  Curve  for  Cantilever  Plate 


Therefore, 


impm  wt  » 


V 

c 


/m 


c 


i 


where  the  constant  is  determined  from  calculations  or  experiment  as  dis- 
cussed above. 

If  this  equation  is  plotted  in  v vs.  M coordinates,  Figure  13, 
the  curve  will  divide  the  plane  into  safe  (below  the  curve)  and  unsafe  (above) 
regions.  For  a particular  impact  situation,  the  values  of  v and  M may  be 
computed.  If  the  point  on  the  graph  corresponding  to  these  values  falls  be- 
low the  curve,  structural  failure  due  to  bending  will  not  occur;  if  the  point 
falls  above  the  curve,  failure  will  occur. 

For  example,  it  is  desired  to  determine  whether  a 0.54-kg  bird  (M  = 1.2) 
impacting  the  blade  at  25  m/sec  will  cause  failure.  Computing  the  dimension- 
less velocity,  we  obtain 


v = (25) 


(0.0071) 
(14.4) (0 . 008) 


= 1.54 


In  Figure  12,  this  is  plotted  as  point  A,  lying  above  the  design  curve  for 
the  cantilever  beam  and  below  the  curve  for  the  plate  based  on  static 
measurements;  thus,  the  former  curve  predicts  failure  whereas  the  latter  does 
not. 
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SECTION  III 


LOCAL  RESPONSE  BY  A FINITE-DIFFERENCE  METHOD 


1.  The  ANEL  Code 

The  ANEL  code  is  essentially  a modified  version  of  the  HEMP  code,  with 
the  addition  of  anisotropic  elastic  constitutive  relation,  and  plane  stress 
problem  capability. 

The  HEMP  code  was  originally  developed  by  Mark  L.  Wilkins  [4]  of  Lawrence 
Radiation  Laboratory,  Livermore,  California.  The  code  has  been  used  success- 
fully to  handle  problems  involving  elastic,  plastic,  and/or  hydrodynamic  constitu- 
tive equations,  In  [4],  example  calculations  included  the  problem  of  the  detona- 
tion of  an  explosive  charge  against  a copper  plate  and  the  detonation  of  an 
explosive  charge  in  an  iron  cylinder.  Karpp[10]has  also  used  HEMP  success- 
fully for  such  classical  problems  as  the  expansion  of  an  idea  gas  under  high 
pressure  in  a tube,  the  impact  of  two  identical  blocks,  and  the  problem  of  a 
cylindrical  explosive  charge  accelerating  an  aluminum  disc. 

HEMP  is  described  as  a two  dimensional  finite  difference  "Lagrangian" 
code  since  it  traces  the  position  of  the  material  particle.  In  the  elastic 


range  the  material  is  assumed  to  follow  the  isotropic  Hooke's  lawand  in  the 
plastic  range  it  is  assumed  to  obey  von  Mises'  yield  condition.  Two  general 


classes  of  problems  can  be  solved  by  the  code,  the  plane  strain  problem  and 
the  axisymmetrical  type  of  problem  where  the  x axis  is  taken  as  the  axis  of 


symmetry.  The  code  itself  is  well  documented  by  E.D.  Giroux  [11]. 

The  modification  of  HEMP  involved  the  replacement  of  the  isotropic  Hooke's 


law  by  a generalized  linear  anisotropic  constitututive  equation.  The  equation 
of  state  was  replaced  by  the  anisotropic  elastic  relation 


°m  ” 3 [Ciiki  eki  + Ciikk  Em] 

which  is  an  equation  relating  the  spherical  component  of  stress,  o^,  to  the 
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mean  strain,  e . In  order  to  give  HEMP  the  ability  to  handle  anisotropic 
m 

plane  stress  problems,  some  changes  to  the  strain  displacement  relations 
had  to  be  effected.  With  these  modifications,  the  new  code  ANEL  can  handle 
any  type  of  anisotropic  linearly  elastic  material  and  hence  any  type  of  com- 
posite layup  for  the  plane  strain  and  plane  stress  cases.  Because  of  the  in- 
herent nature  of  the  axisymmetrical  type  of  problem,  this  aspect  of  the  code 
is  limited  to  material  symmetries  of  no  greater  complexity  than  transversely 
isotropic. 

In  the  next  section  the  governing  equations  used  in  ANEL  will  be  presented. 
Also  a brief  description  of  some  of  the  finite  difference  techniques  will  be 
given.  .The  remainder  of  the  chapter  will  deal  with  the  various  test  and  sample 
calculations  made  with  ANEL. 
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2.  Governing  Equations  and  Numerical  Scheme 

The  ANEL  code  is  basically  the  HEMP  code  with  modified  anisotropic 
elastic  constitutive  equations.  In  addition,  we  have  also  added  in  the  ANEL 
code  the  capability  of  treating  plane  stress  problems. 

The  formulation  of  the  equations  and  the  numerical  scheme  used  in  the 
HEMP  code  have  been  presented  by  Wilkins  [4]  and  Karpp  [10].  Wilkins'  code  and 
presentation  are  for  the  one-dimensional  problem  (uniaxial  strain),  the  plane 
strain  problem,  and  the  axisymmetrical  problem.  Karpp' s presentation  of  the 
governing  equations  is  in  Cartesian  coordinates  only.  For  clarity  and  con- 
tinuity, we  shall  present  all  the  governing  equations  and  numerical  scheme  of 
ANEL  independently;  although  many  of  the  equations  are  identical  to  HEMP. 

All  the  governing  equations  used  are  written  with  respect  to  a coordinate 
system  fixed  in  space,  which  will  be  called  the  space  coordinates.  All  vector 
and  tensor  quantities  in  the  governing  equations  are  written  with  components 
in  these  space  coordinates.  ANEL  may  be  called  a Lagrangian  code  because 
the  properties  and  positions  of  the  mass  particles  are  traced  in  time.  However, 
the  position,  velocity,  stress  and  strain  of  these  particles  are  all  written 
with  respect  to  the  space  coordinates.  in  this  sense  the  code  is  Eulerian  in 
nature.  A better  description  would  be  to  call  the  code  a hybrid  Eulerian- 
Lagrangian  code. 

The  conservation  of  momentum  equations  used  are 


P 


P 


D2x, 


Dt 

D2x„ 


Dt 


0u.i  + ’,(2-T) 


12 


(^22  ~ °33^ 

CT2j , j + ~~ 


(41) 


(42) 
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Here  the  index  j takes  on  the  values  1,2;  comma  j represents  differentiation 
with  respect  to  the  coordinate  x ^ . The  material  derivative  taken  along  a 
particle  path  is  represented  by  D/Dt.  For  a plane  stress  or  plane  strain 
problem  the  space  coordinates,  x-^»x2  would  be  the  Cartesian  coordinates  x,y. 

For  an  axisymmetrical  problem  the  space  coordinates  x^.x^  would  be  the  cylin- 
drical coordinates  z and  r.  The  position  of  the  particle  is  given  as  x^jX^. 

The  stress  tensor  is  represented  by  cr„  and  the  density  by  p.  The  tracer  y 
has  the  following  definitions:  y = 0,  plane  strain  problem;  y = 1,  axisymmetrical 
problem;  and  y = 2,  plane  stress  problem. 

The  conservation  of  mass  equation  is  given  as 


m + 0 ;i,i + py<2-y)  -r2  * ‘ A *1,1 

+ B itizil  B ;2>2  + P xtatt  c(i2>1  ♦ i1>2, 


(43) 


Here  the  index  i takes  values  of  1,2;  x^  and  x^  are  the  velocities  of  the 
particle  in  the  x^  and  x^  directions,  respectively.  The  last  three  terms  are 
needed  only  for  the  plane  stress  case  and  are  discussed  in  Appendix  C . 

The  strain  rate  e is  related  to  the  particle  velocity  by  the  equations: 


i 


'll  ■ 2<*i,J  + *j,l> 


i,j  = 1,2 

= y(2-y)  — + TlT-1?  A 


'33 


xillil  B * + T(y-l)  +x  ) 

2 B 2,2  4 U2,l  1,2; 


(44) 


(45) 


• = y.Cy.-ii  D • + ylir-iJ.  E x 

e13  2 1,1  2 2,2 


+ 1-(X-~1)  + * ) 

4 VX2,1  xl,2; 


(46) 
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• = Illzil  G x + iGCil  R X 

23  2 1,1  2 2,2 


+ Illzil  T(*  + * ) 

4 UX2,1  Xl,2; 

Once  again  the  coefficients  A,B,C,D,E,F,G,R,  and  T are  needed  only  for  the 


plane  stress  case  and  will  be  discussed  in  Appendix  C. 

The  strain  rate  deviator,  e is,  as  usual,  defined  by 


e.  6 ..  e 

ij  ij  iJ  m 


i.J  = 1.2,3 


where  the  scalar  e is  given  by 
m 


e = — e 
m 3 qq 


q = i.2,3 


and  6^  is  the  "Kronekar  Delta"  (6^  = 1 if  i = j and  = 0 if  i 4 j). 
The  stress  deviator  S is 


Sij  °ij  6ij  am 


i,J  = 1,2,3 


where  the  scalar  a is  given  by 
m 


a = — a 

m 3 qq 


q = 1,2,3 


The  constitutive  relationships  are  written  in  terms  of  the  stress 


deviator,  S;  the  spherical  component  of  stress,  o , the  strain  rate  deviator, 
e;  the  mean  strain  rate,  e^;  and  the  stiffness  constants, 


= 3 t3CijkS,  “ Cqqkl  6ij]ek*. 


+ ± [ 3C  ,,,  - C , . 6 ] e 

3 ijkk  qqkk  ij  m 


u 2 

~Dt  = 3 [Ciikl  ekf.  + Ciikl  SJ 


where  Is  the  general  stiffness  matrix  of  the  anisotropic  elastic  medium. 
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The  derivation  of  these  two  equations  from  the  generalized  Hooke's  law 


°ij  ^ijkJl  £kJl 


is  given  in  Appendix  D.  Equations  52  and  53  represent  three-dimensional 


equations,  and  the  indices  run  over  1,2,  and  3.  The  time  derivative  oper- 
ating on  the  stress  deviator  tensor  in  Eq. (52)  is  the  Jaumann  [12]  derivative. 


This  derivative  gives  not  only  the  variation  of  a tensor  with  respect  to  time, 
but  also  the  variation  of  the  tensor  due  to  the  rotation  of  the  coordinate 


system  to  which  it  is  referred. 

Equations  (41)- (53)  constitute  29  independent  equations  ((52)  and  (53) 

are  equivalent  to  6 independent  equations),  in  the  29  variables:  a , S„ , 

• • • 

o , £..,  e..,  e , p,  x,,x..  Now  a short  description  of  the  finite  difference 
m ij  ij  m 1 2 

scheme  used  in  ANEL  will  be  presented. 

The  governing  equations  are  integrated  by  a finite  difference  technique. 
The  material  is  initially  divided  into  grid  zones  by  a rectangular  mesh.  The 
sides  of  the  grid  zones  are  initially  parallel  to  the  fixed  space  coordinate 
system  x^,  x^.  Each  grid  zone  always  encompasses  a fixed  amount  of  mass  and 

will  move  with  this  mass  as  it  translates,  distorts,  dilates,  and  rotates.  In 

• • 

other  words  the  grid  zone  is  fixed  in  the  material.  Tensors  such  as  e,  e,  a,S 

and  scalars  such  as  p,  c , and  a are  calculated  as  the  average  over  the  grid 

m m 

zone  area.  Vectors  such  as  x^,  x?,  the  particle  position  and  particle  velocity 
x^,  x2  are  referred  to  the  grid  zone  corners. 

The  material  derivative  D/Dt  is  approximated  in  finite  difference  form  as 


_d_  f ^ fo(t+At)  ~ fo(t) 

Dt  f0  At 


where  fn  is  either  a scalar  or  a vector  fixed  to  the  moving  particle.  The 


y - 


s 


difference  in  f^  is  taken  between  two  successive  times  that  are  separated 
by  the  time  increment  At. 

The  Jaumann  derivative  which  is  used  in  Equation  (52')  operates  on  the  stress 
tensor.  In  finite  difference  form  it  is  approximated  by 


il 


jot 


it  [sij  (t+At)  - Vt)] 


(55) 


The  stress  states  in  Equation  (55)  are  written  with  their  components  in 
the  space  coordinates.  The  stress  components  calculated  at  time  t 

have  rotated  an  angle  during  the  time  At,  and  this  shall  be  designated 

as  S^.(t).  In  Equation  (55)  S_(t)  are  the  projections  of  S!^  (t)  on  the 
space  coordinates,  or 


aik  V SM(t) 


(56) 


where 


ij 


cos(-w21) 

- sin(-a>21 

sin(-w21) 

COS  (-0)  ^ 

0 

0 

0 

0 

1 


(57) 


and 


. • At 

“21  = At  “21  = T 


3X2 

3x, 


(58) 


For  small  deformation,  the  difference  between  (t)  and  S„  (t)  is  negligible. 

The  partial  space  derivatives  are  obtained  by  the  use  of  Green's  Theorem 
and  the  Mean  Value  Theorem.  For  example  to  find  -g^—  we  use  Green's  Theorem 
which  states  that 


3f 

9x 


- dx  dx  = U 
1 J ( 


f dx„ 


(59) 
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in  conjunction  with  the  Mean  Value  Theorem  stating 


As  mentioned  before  when  describing  the  constitutive  equation  52 

the  stiffness  constants  must  be  written  with  components  in  the  space 

coordinates  x, , x0 . The  material  stiffness  constants  C'  , are  referred 
1 ijki 

to  a material  coordinate  system  X^,X2  iixed  in  the  material. 


-*•  x. 


Figure  15.  Material  Coordinates 


Since  an  originally  rectangular  material  grid  may  not  stay  rectangular,  the 
"material  coordinates"  are  not  always  automatically  defined.  We  have 
arbitrarily  assigned  the  X2  axis  to  be  fixed  with  one  side  of  the  material 
grid  connecting  points  (J,k)  and  (J,k+1);  as  a result,  the  X^  axis  which 
is  orthogonal  to  X,,  is  not  necessarily  along  the  side  of  the  grid  connecting 
(J,k)  and  (J+l,k)  ( See  Figure  2).  The  angle  between  the  X2  axis  and 

the  vertical  is  given  by 


0 = tan 


-1  J xi(J»k+1)  “ x1(J,k)| 
x2(J,k+l)  - x2(J,k)| 


(65) 


The  stiffness  constants  can  be  written  with  components  in  the  space  coordinate 
system  by  the  transformation 


Ci  jkS. 


A . A . A A C 
mi  pj  qk  r£  mpqr 


where  the  direction  cosines  A 


ij 


are  given  by 


A 


ij 


COS0 

sin0 

0 


- sin0  J 

cos0  0 

0 1 
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(66) 


(67) 
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In  the  HEMP  code  only  the  strain  rate  tensor  is  calculated;  strain  is 
never  defined.  The  strain  increment  de^  is  referred  to  the  current  dimension 
and  can  be  called  the  natural  strain  increment.  In  ANEL  the  strain  has  been 
calculated  by  the  use  of  the  following  finite  difference  approximation. 


:..(t+At)  = e../t+At^  At  + e . . 

iJV  T 7 


Here  the  mass  particle  has  been  followed  and  strain  has  been  obtained  by  inte- 
grating step-by-step  in  time.  This  may  be  considered  as  a material  time 
integral  as  defined  by  Malvern  [13]  (p.  151).  This  strian,  however,  is  clearly 

an  approximation  because  no  contribution  due  to  the  rotation  of  the  mass 
particle  as  it  advances  in  time  has  been  considered.  The  strain  rate  is 
always  referred  t o the  space  coordinate  system.  A more  correct  way  to  get 
the  strain  would  be  to  treat  it  the  same  way  as  was  done  with  the  stress  tensor 
by  taking  into  account  the  rotation  of  the  mass  particle.  In  =hort  the  values 
of  strain  that  have  been  calculated  in  ANEL  may  be  considered  as  Lagrangian 
with  respect  to  translation  but  Eulerian  with  respect  to  rotation. 

The  set  of  algebraic  finite-difference  equations  are  used  to  calculate 
from  the  known  quantities  at  time  t the  solution  of  tne  unknowns  at  time  t+At. 
Figure  16  gives  a flow  chart  of  the  order  in  which  the  governing  equations 
(in  finite  difference  form)  are  used.  Note  that  we  begin  with  the  stress 

state  a..,  tne  density  p,  and  tne  particle  position  x, ,x„  known  at  t = t. 
ij  1 2 

Eqs- (41)-(42)  are  used  to  obtain  the  particle  velocity  half  way  between  the 

two  time  planes  or  t = t + Here  we  have  really  integrated  Eqs. (41)- (42) 

Dx,  L 

i • • 

once  in  time  to  get  (also  written  as  x^,  -x.^  ) . By  using  the  following 

finite  difference  equation  the  particle  position  at  t = t + At  is  obtained 


x. (t+At)  = x. (t+  At)At  + x. (t) 
1 3-  “2 
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at  t = t + At 
from  Eq . (50) 


Artificial 
Viscousity  ? 
Eq.  (69) 


Figure  16.  Flow  Chart  of  Calculations  in  ANEL  (Cont.) 


are  obtained.  Next 


S 


Now  using  Eqs  • (65)-(67  ) the  stiffness  constants 
the  type  of  problem  considered  is  decided:  plane  strain,  y = 0;  axisymmetric, 
y = 1;  and  plane  stress,  y = 2.  Note  only  for  the  case  y = 2,  will  the  coeffi- 
cients A,B,C,D,E,F,G,R,  and  T be  calculated  as  described  in  Appendix  C.  The 
next  group  of  steps  produce  the  density  p at  t+At  by  use  of  Equation  (43),  the 
strain  rate  tensor  at  t+^p  from  Eqs. (44)-(47) , and  the  angle  of  rotation 
^21  at  t+At  from  Equation  (58).  Now  knowing  , Eqs. (48)-(49)  are  used  to 


At 


give  e and  e at  time  t + . 
ij  m 2 


As  is  indicated  in  Figure  16  the  strain 


tensor  e at  time  t+At  is  calculated  by  Eqs.  (68). 


S . . and  o are  obtained 
ij  m 


at  time  t+At  by  using  the  constitutive  equations  Eqs.  (52)-(53).  Equation  (50) 


now  gives  at  time  t+At.  The  artificial  viscousity  is  calculated  as 
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where 


1 

v = — 


P 

• • 

v = _ £_ 

V P 

C = a constant 
A = grid  zone  area. 

This  term  is  subtracted  from  a.,  and  and  a . The  purpose  of  this  term 

11  zz  m 

is  to  smooth  out  the  effect  of  a shock  wave  over  several  grid  zones  and  is 
based  upon  the  von  Neumann-Richtmyer  method.  The  time  is  now  incremented 
by  At  and  if  this  new  time  is  not  greater  than  the  stopping  time,  tgt  the 
program  repeats  the  process.  It  should  be  noted  that  the  time  increment  At 
is  an  inputted  constant  in  ANEL. 

In  the  next  section  several  runs  and  comparisons  using  ANEL  are  discussed. 
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3.  An  Axisymmetrical  Impact  Problem 

As  a partial  check  of  the  code,  a problem  with  isotropic  material  was 

solved  by  ANEL  and  the  results  compared  with  those  of  HEMP.  The  problem 

involved  au  isotropic  cylinder  subjected  to  a point  load  along  its  axis  of 

symmetry.  The  cylinder  had  a length  of  4 cm  and  a radius  of  4 cm.  Young's 

modulus  used  was  E = 2.07  Mbar  and  the  shear  modulus  was  G = .796  Mbar.  The 

3 8 

density  of  the  body  was  p = 7.8  gm/cm  . The  point  load  of  2.62x10  dynes 
was  applied  gradually  in  time  and  then  held  constant.  Furthermore,  in  order 
to  avoid  any  instability  in  either  program  due  to  the  discontinuous  nature 
of  a point  load,  the  force  was  replaced  by  a conical  pressure  distribution. 

The  form  of  this  pressure  distribution  in  time  t and  radius  r is 

- ■ -oC1  - a ft  -(t)2) 

p-po[1-fJ 

P = 0 

where  a = 0.5  cm,  tQ  = 1 ysec,  and  Pq  = 0.001  Mbar* 

Both  programs  were  run  as  axisymmetrical  cases.  The  mesh  size  used  in 

the  region  closest  to  the  axis  of  symmetry  was  1x3  mm  and  the  constant  time 

step  used  throughout  the  calculations  was  At  = 0.02  ysec. 

Figure  17,  shows  a time  history  of  the  normal  stress  a at  the 

zz 

location  r = 1.02  cm,  z = 2.12  cm.  The  curve  plotted  is  the  same  for  both 
the  HEMP  run  and  the  ANEL  run.  This  indicates  that  ANEL  can  handle  an  isotropic 
elastic  problem  and  gives  very  favorable  results  when  compared  with  HEMP. 

Another  problem  was  solved  by  ANEL  and  the  results  compared  with  an 
existing  numerical  solution.  The  problem  consisted  of  a semi-infinite  compos- 
ite subjected  to  an  impact  of  a rigid  sphere.  In  Ref.  [14],  Greszczuk  solved 

* ANEL  uses  the  c.g.  ysec  system  of  units:  length  - cm,  mass  - g,  time  - ysec, 
force  - 10l2  dynes,  and  pressure  - Mbar. 


t < t0 
r < a 

t > t0 

r < a 
r > a 
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this  impact  problem  by  a quasi-static  approach.  The  contact  pressure  dis- 
tribution was  calculated  by  the  Hertzian  method,  but  the  response  in  the 
composite  is  calculated  only  as  a static  problem,  inertia  force  was  neglected. 
He  assumed  a surface  indentation  and  obtained  the  pressure  distribution  of  the 
form 


P = q0  t1  " a 

P = 0 r >_  0 

where  a is  the  radius  of  the  area  of  contact  between  the  sphere  and  the  target 
surface  and  qQ  is  the  maximum  pressure  occurring  at  the  center  of  this  area 
of  contact.  This  approach  is  shown  schematically  in  Fig.  18.  By  the  use  of 
the  finite  element  program  SAAS  III  he  was  able  to  get  the  local  stress  dis- 
tribution in  the  composite  target  due  to  the  above  pressure  distribution. 

For  the  comparison  between  ANEL  and  Greszczuk's  work  the  composite  chosen 
was  a transversely  isotropic  glass  epoxy  laminate.  The  elastic  constants 


used  in  c.g.  psec  system  of  units  are 


E 
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= 0.2  Mbar 
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=0.115  Mbar 

was 

p = 1.995  gm/cm3. 

The 

stiffness  constants 

C. ..  „ as  calculated  from  the  above  elastic  constants  are  given  in  Table  7. 
ij  ki- 
ln [14]  the  maximum  pressure  and  radius  of  contact  were 


q^  = 0.00031  Mbar 
a = 0.5  cm 

In  the  ANEL  calculation  the  pressure  distribution  was  applied  gradually  in 
time  and  then  held  constant  as  shown  in  Fig.  19.  A true  Hertzian  pressure 
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Sphere  on  a Semi-Infinite  Anisotropic  Medium 
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Table  7.  STIFFNESS  CONSTANTS  FOR  THE 
COMPOSITE  USED  IN  THE  HERTZIAN  IMPACT  PROBLEM 
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distribution  in  time  is  included  on  this  figure  for  comparison.  The 
loading  pressure  distribution  used  in  ANEL  is  now  summarized  below  as 
P = qf 


*o  ll1  - f t 

P = 0 


r < a 
t < t 


0 


r < a 

c > tQ 


r > a 


where  t^  = 1 psec 

The  target  was  a cylinder  of  radius  4.24  cm  and  depth  4.24  cm.  The  grid 
zoning  used  is  shown  in  Figure  19.  The  constant  time  step  used  was 
At  = 0.02  usee. 

Figure  20  gives  a time  history  of  the  normal  stress  o^/q^  at  a 

location  on  the  axis  of  symmetry  0.9  cm  away  from  the  center  of  impact. 

The  maximum  value  of  normal  stress  occurs  at  t = 4.4  usee. 

Figure  21  is  a plot  of  the  axial  distribution  of  normal  stress.  The 

solid  line  is  from  ANEL  and  was  plotted  for  the  time  t = 4.4  usee.  The  dashed 

line  is  the  static,  finite  element  solution  from  [14].  Although  both  solutions 

start  at  the  same  value  of  stress,  o /q_.  = 1 at  z/a  = 0,  the  dynamic  solution 

z 0 

is  higher  at  all  other  points  than  the  static  solution. 

Figure  22  is  a plot  of  the  radial  distribution  of  shear  stress  along 
the  radial  line  z/a  = 0.6.  The  solid  curve  is  plotted  from  the  ANEL  solution 
at  t = 2.4  usee.  The  maximum  value  of  shear  stress  along  this  radial  line 
occurs  at  this  time.  The  dashed  line  is  the  static  finite  element  solution 
of  [14].  The  peak  shear  stress  from  both  approaches  occurs  approximately  at 
r/a  = 0.8.  However,  the  peak  from  the  dynamic  solution  is  approximately  fifty 
percent  higher  than  the  static  solution. 
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Figure  20.  Time  History  of  Normal  Stress  at 


Figure  22. 


Radial  Distribution  of  Shear  Stress  at  2. A usee  and  z/a=0.6 


Glass  Epoxy  Under  Impact  of  Rigid  Sphere 
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4.  A Plate  Edge  Impact  Problem 

Another  calculation  was  made  of  a cantilever  plate  impacted  by  an 
aluminum  projectile  on  its  edge.  Experiments  for  this  impact  case  were  also 
conducted  and  the  results  from  the  two  were  compared.  Figure  23  shows 
the  experimental  set  up.  The  plate  was  7 in.  (17.78  cm)  wide,  13  in.  (33.02  cm) 
long,  and  0.28  in.  (0.71  cm)  thick  and  was  rigidly  supported  at  its  upper  edge. 
The  plate  was  made  of  50  layers  of  graphite  epoxy  and  the  lay-up  of  the  laminae 
were  ( (90, 45 ,90, -45) ^ , 90) ^ where  the  angle  was  measured  with  respect  to  the 
x-axis  which  was  defined  as  the  lower  edge  of  the  plate.  The  material  con- 
stants for  each  lamina  were 

E11  = 17>7  x 1C)6  PSI  v23  = 0,2 

E22  = e33  = !-3  x 1q6  PSI  g12  = Gi3  =0’55  x 1q6  PSI 

V12  = V13  = 0,3  G23  =0,54  X 106  PSI 

where  the  subscript  1 was  the  fiber  direction  and  2,3  were  the  directions  per- 

3 

pendicular  to  the  fibers.  The  density  of  the  plate  was  p = 1.51  g/cm  . The 
plate  was  impacted  by  an  aluminum  projectile  of  mass  0.06  lb  (27  g)  which 
travelled  at  a constant  velocity  of  40  ft/sec  (12.2  m/sec).  The  path  of  the 
projectile  was  parallel  to  the  x axis.  The  point  of  impact  was  on  the  left 
hand  edge  at  a distance  of  2.5  in  (6.35  cm)  above  the  x axis.  At  the 
coordinates  x = 5 in.  (12.70  cm),  y = 2.5  in.  (6.35  cm)  two  strain  gages  were 
placed,  one  on  the  upper  side  and  one  on  the  lower  side  of  the  plate,  to 
monitor  bending  stress. 

The  conventional  method  of  handling  a composite  plate  problem  is  to 
use  the  lamination  theory  (see  ref.  [9]).  The  known  values  of  the  lamina 
engineering  constants  E’s  and  v'3  are  first  converted  into  plane  stress 
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"reduced"  stiffness  constants  Q...  These  Q. . are  then  transformed  into  the 

13  iJ 

plate  coordinates  and  summed  to  give  the  A„ , and  D_  matrices.  For  the 

present  case,  the  lay-up  is  symmetric  and  bending  is  not  induced;  therefore 


only  A. . are  needed. 
iJ 

Due  to  the  arrangements  in  HEMP,  we  found  it  more  convenient  to  express 
the  stress-strain  relations  in  a different  way  for  ANEL.  The  layered  com- 
posite is  first  "smeared"  into  an  equivalent  homogeneous  anisotropic  material 
with  a set  of  engineering  constants  and  a corresponding  set  of  stiffness  con- 
stants (contracted  notation),  both  three-dimensional.  These  three- 

dimensional  are  to  be  fed  into  ANEL.  The  plane  stress  assumptions  and 
equations  are  already  incorporated  in  ANEL. 

The  plate  constants  A^  can  also  be  obtained  from  the  three  dimensional 
of  the  smeared  material.  Substituting  these  C„  into 


0 , c - Cl3  CJ3 

13  13  C33 


i,j  = 1,2,6 


yields  the  reduced  stiffness  Q for  the  composite,  which  after  multiplying 
by  the  total  thickness  of  the  laminate,  gives  exactly  the  A as  obtained 
from  Eq.  (3-25)  of  Ref.  [9]. 

For  the  present  laminated  plate,  the  "smeared"  three-dimensional 
engineering  constants  are 
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The  stiffness  constants  in  Mbars  for  this  equivalent  medium  are  given  in 
Table  8.  The  material  symmetry  of  the  equivalent  medium  is  orthotropic. 

In  the  numerical  calculation,  the  plate  is  divided  into  square  grids, 
each  with  sides  of  0.47  in.  (1.2  cm).  The  lower  edge,  left  hand  side  edge, 
and  the  right  hand  side  edge  were  specified  as  free  boundaries.  The  upper 
edge  was  specified  as  a fixed  boundary. 

In  order  to  model  the  impact  loading  the  Hertzian  method  described  by 
Goldsmith  [15]  was  used.  First  the  projectile  was  simulated  by  an  aluminum 
ball  with  the  same  mass  and  a radius  of  1.3  cm,  and  the  same  velocity  as  the 
original  projectile.  The  pressure  distribution  according  to  Eqs.  (4.13,  4.33) 
of  [15]  was  then 


P = P 
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t < t 
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r < a 
— m 


where 


P = 0 


P = 
m 
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2ir  a 
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and/or  C — C0 

r > a 
— m 


and 


a = the  maximum  radius  of  the  area  of  contact  = 0.21  cm 
m 

tg  = the  impact  duration  = 81  ysec 

12 

F = maximum  impact  force  = .0014x10  dynes  . 
m 

Because  of  the  large  mesh  size  used,  this  pressure  distribution  covers  less 
than  one  grid  zone  and  was  considered  unsatisfactory;  instead,  a rectangular 
pressure  distribution  of  the  same  duration  in  time  with  the  same  total  force 
but  larger  contact  area  was  used,  and  is  given  below  as 


j 
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P = Pq  sin 


c - t0 


4.7cm  8.3  cm 


P = 0 


t > tg  and/or 


4. 7cm  > y > 8.3cm 


-.-4 


where  Pq  = 5.4x10  Mbar.  As  can  be  seen  this  pressure  distribution  spreads 


over  three  grid  zones. 

Figure  24  shows  the  comparison  between  ANEL  and  experiment  of  the 
strain  at  the  location  of  the  strain  gage  at  x = 5 in.  (12.70  cm)  and 
y = 2.5  in.  (6.35  cm).  Great  care  was  taken  in  the  experimental  study  to 
obtain  identical  traces  for  the  gages  located  on  the  upper  and  lower  sides 
of  the  plate.  Any  deviation  between  the  two  would  indicate  bending  and 
hence  the  plane  stress  assumption  would  not  be  applicable.  The  dashed  line 
in  the  figure  represents  the  reading  of  strain  from  the  strain  gage.  About 
ten  experimental  trials  were  run  and  a statistical  study  was  done  on  the 
peak  strains  and  their  times  of  occurrence  and  also  on  the  time  that  the 
strain  curve  crossed  the  time  axis.  The  coefficient  of  variation  (i.e.  the 
standard  deviation  divided  by  the  mean)  for  these  readings  was  found  to  be 
no  greater  than  12%  in  all  cases.  This  indicates  good  reproducibility  of 
the  experimental  results.  The  solid  curve  represents  the  solution  obtained 
by  ANEL.  As  can  be  seen  the  analytical  solution  has  the  same  shape  as  the 
experimental  results  and  the  times  of  the  peaks  and  zero  strain  are  in  good 
agreement.  However,  the  magnitude  of  the  peak  strains  from  ANEL  are  greater 
than  those  obtained  from  the  strain  gage  readings.  This  difference  is  about 
110%  for  the  first  peak  occurring  at  65  psec  and  about  130%  for  the  second 
peak  occuring  at  120  psec. 

There  are  several  explanations  for  this  large  discrepancy  in  peak  strain 
shown  in  Figure  24.  First  the  loading  function  could  be  inaccurate  be- 

cause the  projectile  was  modelled  as  a ball  of  radius  1.3  cm  when  in  reality 
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Figure  24.  Comparison  for  Edge  Impact  of  Composite  Plate 
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it  was  a cylinder  2.54  cm  long  and  2.16  cm  in  diameter.  The  striking  end  of 
this  projectile  was  rounded  with  a radius  of  2.54  cm.  Perhaps  by  using  either 
the  radius  of  the  rounded  end  for  the  radius  of  the  ball  or  the  diameter  of 
the  cylinder  for  the  diameter  of  the  area  of  contact  would  have  given  a more 
realistic  loading  function.  Another  reason  for  the  discrepancy  between  the 
computer  study  and  experiment  is  that  ANEL  is  a completely  elastic  code.  No 
provision  was  made  for  the  plastic  deformation  during  impact.  Likewise  the 
energy  dissipation  and  wave  dispersion  were  not  accounted  for  in  ANEL.  The 
wave  propagates  from  the  impact  point  to  the  strain  gage  in  a "weak" 
direction  of  the  laminated  plate,  i.e.  it  is  90°  or  45°  from  the  fiber  direc- 
tion; the  wave  is  essentially  travelling  thru  the  matrix.  It  can  be  expected 
that  the  viscoelastic  effect  was  strong  and  much  of  the  energy  was  dissipated 
before  reaching  the  gage. 

Figure  25  shows  a plot  of  the  strain  across  the  plate  from  the  left 
hand  edge  to  the  right  hand  edge  at  y = 1.87  in.  (4.75  cm)  and  t = 48  sec. 

The  strain  that  occurs  at  the  left  hand  edge  on  this  curve  is  the  maximum 
value  of  strain  occurring  in  the  140  usee  for  which  ANEL  was  run.  As  seen 
the  strain  decreases  as  one  moves  across  the  plate  from  left  to  right. 

In  summary  one  can  say  that  ANEL  tends  to  give  higher  readings  than  ex- 
periment but  does  model  the  shape  of  the  wave  fairly  accurately.  More  work 
must  be  done  with  the  code  in  order  to  introduce  plasticity.  Also  some  addi- 
tional study  must  be  performed  because  in  the  plane  stress  case  of  ANEL  some 
numerical  anomaly  appears  at  the  free  boundary.  However,  as  ANEL  exists  now, 
one  can  say  that  the  program  does  at  least  give  conservative  results. 


Strain  Distribution  Across  Composite  Plate  at  y = 1.87  in.  and  t = 48  usee 
Maximum  Strain  at  x = 0.47  in.,  y = 1.87  in.,  t = 48  usee. 


APPENDIX  A 


STATIC  SOLUTION  OF  A SIMPLY  SUPPORTED  ORTHOTROPIC  PLATE 


In  this  appendix,  we  will  present  the  details  of  the  derivation  of  the  con- 


stants K1  and  d12  for  the  central  transverse  impact  of  a simply  supported 
orthotropic  plate.  As  explained  in  the  main  text,  these  constants  are  needed 
for  developing  the  design  curve. 


The  static  deflection  of  such  a plate  due  to  an  arbitrary  load  distribution 
p(x,y)  is  [Al], 


w(x,y)  = l l — 

ra=l  n=l  m tt 
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mn a b 
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4 I , v nrnx  . niry  , 

p = — p(x,y)  sin  — sin  -rf-  dxdy 

mn  ab  JQ  JQ  a b 


For  a centrally  applied  point  load  P, 
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The  maximum  deflection  occuring  at  the  center  of  the  plate  is. 
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where 


we  have 


The  tensile  strain  in  the  x-direction  due  to  bending  is 


which  will  have  its  maximum  value  at 


We  will  derive  the  relationship  between  the  strain  at  this  point  e.  and  the 


maximum  deflection 


where 


Substituting  for  P,  we  obtain 


so  that 


The  approximation  assumption  used  above,  Eq.  (A2) , has  been  suggested  by 
Timoshenko  and  Woinowsky-Krieger  [Al]  for  simplifying  the  mathematics  of  several 
plate-bending  problems  and  is  exact  for  the  case  of  an  isotropic  plate.  The 
advantage  of  this  assumption  is  to  reduce  the  parameters  for  describing  the 
geometry  and  anisotropy  of  the  plate  to  a single  quantity,  n.  For  the  composite 
plate  on  which  the  impact  experiments  were  performed,  the  computed  values  are 

H = D, „ + 2D  = 13140  lb-in. 
oo 

and  /D  D = /d  D = 11517  lb-in. 
x y _L-L  Z.2. 

a difference  of  14%. 


Reference 

Al.  Timoshenko,  S.,  and  Woinowsky-Krieger,  S.,  Theory  of  Plates  and 
Shells,  Second  Edition,  McGraw-Hill  Co.,  New  York,  1959 
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APPENDIX  B 

COMPUTER  PROGRAM  FOR  CALCULATING  TIMOSHENKO  SOLUTION 
OF  PLATE  IMPACT 


This  appendix  presents  a listing  of  the  FORTRAN  computer  program, 
named  SMINC3,  for  calculating  the  Timoshenko  small- increment  solution  of 
the  central  transverse  impact  of  a rectangular  orthotropic  plate,  Eq.  (29). 

The  program  as  listed  is  set  up  to  solve  an  impact  problem  having  the 
following  parameters: 

v = 2.45  m/sec 

axbxh  = 171  x 349  x 7.1  mm 


m^  = 0.660  kg 

m2  = 0.816  kg 


The  plate  bending  stiffnesses,  computed  according  to  the  Whitney-Pagano 

anisotropic  plate  theory  [9],  are: 

D = 680.4  N-m 

Dl2  = 473.2  N-m 

D22  = 2488.8  N-m 

D,,  = 505.7  N-m 
66 

D,,  = D„,  = - 49.1  N-m 
lb  z0 

A..  = Acc  = 2 . 687xl07  N/m 
44  55 


The  plate  is  treated  as  orthotropic  and  the  and  terms  are  neglected. 

The  program  itself  calculates  the  value  of  the  Hertzian  contact  stiff- 
ness constant  based  on  supplied  values  of  isotropic  elastic  properties  of 
the  plate  and  impactor.  The  values  used  are: 


plate  (Gr/Ep) 


impactor  (steel)< 


E = 9 . 8xl09  N/m2 
v = 0.3 

E = 2.1x10**  N/m2 


[ \>  = 0.33 

The  contact  radius  of  the  impactor  is  25.4  mm. 
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C SMINC3  --  PLATE  IMPACT 

CQMPI  TATIOim  OF  1 ImOqHENkO-TYPE  SmALL“INCREhEnT  solution  OP  hERT/IAN 

c impact  uf  a simply  supported  orthotrqpic  rectangular  plate 

C MODIFIED  SUN  & CHATTOPADHY AY  METHOD  

C _ D T - TIME  INCREMENT 

_C DDT  " N U « Oh  TIME  INCREMENTS  CALCULATED 

C TDL  ’ TOLERANCE  OF  ERROR  IN  NONLINEAR  SULuTIQN 

_C  V ” IMPACT  VELOCITY 

~C  X A / YS,  H - DIMENSIONS  de  plate  IN  X*  Y / L directions 

_c Ml/  M2  - HASS  OF  P L A Tt  > IMPACTOR 

C f.  1»  E 2 " ELASTIC  MODULUS  UF  PLATE/  IMPACTOR 

_C P N 1 / P K 2 ~ PU1SSOnmS  KATIU  UF  PLATE/  1MPACTUR 

C NhM/  NHM  - NUMBER  Of  VIBRATION  MODES  CONSIDERED 

_C K 2 ~ HERTZIAN  CONTACT  STIFFNESS 

~C  XNX  / Y N Y “ INITIAL  STPESStS  IN  X/Y  DIRECTIONS 

PEAL  L/H1/H2/K2/M1/L12/L13/L22/L23/1  33 

DIMENSION  F ( 5000  ),  (i  ,<  C 5000  ) /GXX(  5000) /GYY(  5000) 

N H M = 2 5 ” 

NHN=25 


NDT=lUO 
TOL  = 1 .OE-05 
0T  = 1 .0E“06 

X N X = 0 « 

Y N Y = 0 . 

Dll=6O0.4 


012=473.2 
022  = 2486.  t 


PRIM  1 
CALCULATE  M 


D09CK= 
Ga(*  ) = 
G X X ( is  ; 
VO  GYY(K) 

0 0 2 A A 

T £ K P.  A = 

TERM Ac 
L 1 3 = A*: 

DC12A7.' 

T t R M o 5 

T E R M b c 

l 1 1 =U 1 

L 1 2=  C L 

L22  = U< 

L23=A< 

L 3 3 = C A 

C'  = L1  H 
OF  T E hi* 

A C = t L 1 

PC=(Li 

G M E G A i 
Q M E 0 A - 

1 F ( ( M , 
1 U 0 2 F 0 F M A 1 

TERM1 = 
C 2 = 1 .1 
DU?4 (* 

C 1=12 

C 2 = E U i 

A A A = ( ( 

CONSTRUCT  T / 
G R ( IS)  = 

g"  x x ( t\ : 
2A7  cyy(k; 
p=b  + l < 
e = i . /. 
S = 2 • / ; 

0 U 1 u u r 

"CREATE  HEAD] 
. F (Min 
1 0 0 3 F OR l-  m i 
PRINT  1 

print i 
70  T=N*ul 
a = v * r 

I f ( Is  . t 

BUM=0 , 

MM  1 = N" 

S 0 M 1 = ( 
Sum  = g i 

0Cj5j="i 

K = N ' * J 

R = K 

Sum  i =f 

SUM  = G(. 
5 B U M = ti  l 


1001 
a r u « a 


-l,NU 

=JJ  . 

) =0 . 

) = 0. 

= M « P I 
2=1  ER 
bb  ♦ I E 

nl  = 1 * N 
= N » P I 
2=1  EK 
1 1 «TE 
012  + 0 
66  » T E 
4M*T  E 
ASb  + X 
+ L22- 
M = M * L 
1 2*L2 

1 2 * L 1 

2 = 0E  T 
= 6hR  1 
• EU.l 
T ( 11 

= OT*U 
0 

K = 1 > is 


l ^REwufNCits  of  plate 


T 


!*  D66«TERMB2+A55 

> T E R A * T t R M B 

> + U2  2»TERMB2  + AA» 

■ TER  A2  + EA4A  + YNY  )*TERMB? 

i*2 

?.*L12»L23*L13*L22»L13*«2“L11*L23**2 

!2*Ll3)/0*TERMA 

■1*L?3)/U*TEKMB 

'0/P 

■ GA2  ) 

.0.  (n.EQ. 1 ) ) PR  I NT  1002* OMEGA 

. G A * 1 * 1 = "#115,6)  


1 * 1 E R M 1 ) 

•12  )/UHEGA? 

.es  u f summation  functions 

I ( K ) + Am  A 

1 X X C N ) + A A A * A C 

»r  T ( N ) + A AA*B  0 

1(1") 


. » NOT  _ 

iS  F uR  D AT  A AT  T UP 
1*61)  . NE  • 0 ) GO  TO  7 0 
.Ml  ) 

) 3 
11 


E 0 a 1 ) U(J  T U 1 4 


EVERY  PAGE 


1 * N 1 1 


F ( J ) - S U M 1 
U IT  + R * S U M 1 
uM  + t ( J ) * G v.  ( K + 1 ) 


,.r  » 


CHtCK  THAT  F l 1 ) IS  P 0 S>  I T 1 - 

IF  (F  IM£W.GT  « U . JGUTfJl  0 

FNfc.n=“KNLM 

GUTblJ  


FPFns3.U«P*«3*FN*»;j+2.0*B»FN  + 3»u*  A»*2»  P 
6?  F Ntw  = FfJ-F  F N/FPFN 

FRKLiK  = ABSC(FNti.-F'Mj/FiMj 

IF  (LKrtUN.oT.TUL)GuTUl3 

IF ( N . GT . 1 IGuTulO  


10  F (N)=F NEW 

I F ( F N E W . L T . 0 . ) FIN)  = 0» 


CALCULATE  STRAIN  AND  DEFLECTION 
W = 0 . 


PSlXX=U. 

PSI  YY=0. 

PU21 J=1»N 

K=N-J+1  


W=W+F(J)*GW(K) 

PSIXX=PSIXX+FCJ)*GXX(K)  

21  PSI YY=PSIY Y+F C J)*GYYCK) 

W = ri  » C 

FPSx="PSlXX*C*Z 

FPSY=-PSI YY*C*/  


G A M M A = 0 . 

ALPhA  = SIGN( 1 « J » K N F R ) * ( A B S t F~  N E R ) /K2 ) «« S 

PRINT  1000/ N/T,F(N), mI  PHA#h#EPSX,EPSY»faAMMA 

100  ClINT  i in UE ___ 

S TUP 

1 UOO  F0KMAT(1X*I4>7C3X,E15.6))  
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APPENDIX  C 


COEFFICIENTS  FOR  THE  PLANE-STRESS  CASE 


The  terms  containing  the  coefficients  A,B,C,D,E,F,G,R,  and  T in  the 
conservation-of-mass  equation  (43)  and  the  strain-rate  equations  (44)- (47)  are 
different  from  zero  only  for  the  plane-stress  case.  They  result  from  the 
plane-stress  assumptions,  i.e.,  the  stress  components  a.^»  °23’  anc*  °33  eclua'*- 
to  zero  and  there  is  no  variation  in  stress  or  strain  in  the  direction  per- 
pendicular to  the  plane  (i.e.  the  direction).  The  relationship  between 


stress  and  strain  is 


°ij  Cijk«.  £k Z 


Now  since  °23»  an<*  °33  are  zero»  °ne  can  solve  for  e^,  £^3,  and  z^  in 


terms  of  , z^,  and  e ^ and  get 


e33  A E11  + B e22  + C £12 
E13  = D E11  + E e22  + F £12 
£23  = G £11  + R E22  + T E12 


These  coefficients  are  in  terms  of  the  stiffness  constants  j an^ 
Carleone  [16]  have  compiled  these  coefficients  for  various  material  symmetries. 
For  instance  an  orthotropic  material  will  have 


A = - 


C = D=  F = G=  R=  T=  0 


APPENDIX  D 


DERIVATION  OF  CONSTITUTIVE  RELATIONS 


This  appendix  gives  the  derivation  of  the  constitutive  relations  (52) 
and  (53)  which  are  in  terms  of  the  deviator  stress  and  deviator  strain, 
from  the  familiar  generalized  Hooke's  Law  in  terms  of  total  stress  o and 
strain  e.  The  generalized  Hooke's  law  is 

aij  = CijU  £k{,  (D_1) 

where  C. „ are  the  stiffness  constants  of  the  anisotropic  elastic  material. 
ijk£ 

(See  for  instance  [16]  and  [17]) 

The  deviator  stress  S and  the  deviator  strain  e are  given  as  follows: 


S . , = a . - a 6 . . 
ij  iJ  m ij 


e e 6.. 
ij  iJ  m ij 


where 


°m  3 U 

= I 

£m  3 


Note  all  indices  assume  values  1,  2,  and  3. 


Performing  a tensor  contraction  on  Eq.  (D-l)  gives 

°ii  = Ciik{,  ekJl*  (D_6) 

Making  a substitution  of  £ from  Eq.  (D-3)  and  a substitution  of  0^  from 

Eq.  (D-4) , Eq.  (D-6)  then  becomes 


°m  3 | Ciik£  ekR  + Ciikk  em|  (D  7) 

Let  us  now  substitute  Eqs.  (D-2)  and  (D-3)  into  Eq.  (D-l).  After  some 


rearranging  the  following  is  obtained 


S..+o  6..=CJJ1.e,  „+C.  ...  e 

ij  m ij  ijkf.  kS,  ljkk  m 
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Substituting  for  a in  Eq.  (D— 8)  from  Eq.  (D-7)  gives 


" I [3  cy 


- C 


ijk£  qqk£  ijj  k£ 


(D-9) 


+ 3 [3  Cijkk  Cqqkk  6ij]  £m 
Now  taking  the  total  time  derivative  along  a particle  path  of  Eqs.  (D-9) 
and  (D-7)  gives 


^Si1  _ 1 

jy  t 3 

i 

3 


3 C. 


ijk£  qqk£  ijj  k£ 
e 


Do 

i 

Dt 


m 


1 

3 


3 Cij kk  " Cqqkk  6ij 


CiikJl  ekX,  + Ciikk  £m 


m 


(D-10) 


(D-ll) 


where  t is  the  Jaumann  derivative  that  has  been  discussed  in  the  main 


body  of  the  report. 

Eqs.  (D-10)  and  (D-ll)  are  Eqs.  (52)  and  (53)  respectively. 
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LIST  OF  SYMBOLS 

plate  dimension  in  x-direction 

plate  dimension  in  y-direction 

beam  flexural  wave  velocity  = ■/E/p 

plate  flexural  wave  velocity 

plate  flexural  stiffness  matrix 

effective  mass  factor 

flexural  rigidity  of  beam 

contact  force 

depth  of  beam 

radius  of  gyration 

Hertzian  contact  stiffness 

eqivalent  beam  spring  constant 

linearized  contact  stiffness  (spring  constant) 

length  of  beam 

mass  of  structure 

mass  of  impactor 

m /m2 

static  load 
time 

velocity  of  impactor  before  impact 

initial  velocity  immediately  after  impact 

w(x,t)  = beam  deflection 

w(x)  = static  beam  deflection) 

w^(x)  = deflection  of  structure  at  impact  point 

maximum  midspan  deflection 


s 
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